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We present numerical results on electron tunneling in a single-electron box at low temperature. 
The effective action of this device is equivalent to the Hamiltonian of a classical XY spin chain 
with long ranged interactions. Using an efficient cluster algorithm and a new transition matrix 
Monte Carlo approach, we are able to compute the effective charging energy E c in the limit of very 
small tunneling resistance. While previous Monte Carlo simulations were restricted to the weak and 
intermediate tunneling regimes, our method extends the range of _BJ-values by more than 30 orders 
of magnitude. This allows us to clearly observe the exponential suppression of Eq with increasing 
tunneling conductance a. For large, but fixed a, the correction to the leading exponential behavior 
exhibits a crossover from an intermediate temperature behavior at flE c <C 1, to zero temperature 
behavior at (3E C S> 1. We determine this correction in both regimes and compare the numerical 
results to the numerous and controversial theoretical predictions for the strong tunneling limit. 



I. INTRODUCTION 



The so-called single electron box is an example of a dis- 
sipation coupled mesoscopic system with discrete charge 
states. It consists of a low-capacitance metallic island 
connected to an outside lead by a tunnel junction. The 
single electron box has been the subject of numerous 
experimental and theoretica l 2 ' 3 ' 4 ' 5 ' 6 ' 7 investigations, be- 
cause it exhibits Coulomb blockade phenomena due to 
the large charging energy of the island. The presence 
of excess charges influences single electron tunneling and 
this effect is used in single electron transistors 8,9 , which 
can be viewed as a single electron box connected to two 
tunnel junctions for the entrance and exit of electrons. 

The first experimental realization of a single electron 
box was achieved by Lafarge et al} . A circuit diagram 
of this device is shown in Fig. ^ The box with excess 
charge n is controlled by an external voltage source Vq 
to which it is connected through a capacitor Cq and a 
tunnel junction with resistance Rt and capacitance C*. 
An applied gate voltage Vq induces a continuous polar- 
ization charge iiq — CgVg/c The bare charging energy 
Ec = e 2 /2(C t + Cg) sets the energy scale. 
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FIG. 1: Circuit diagram of the single electron box. The box 
with excess charge n is indicated by the dashed rectangle. It 
is connected to a voltage source through a capacitor Cg and 
a tunnel junction with resistance Rt and capacitance Ct- 
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FIG. 2: Illustration of the Coulomb staircase in the limit 
k B T < E c and R t 3> Rk- 
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The classical electrostatic energy for an integer number 
n of additional electrons on the island is given by the set 
of parabolas E n {nc) — Ec(n — iig) 2 , as shown in Fig. [21 
In the ground state, the number of excess charges on the 
island is therefore a staircase function with unit jumps 
at riG = 1/2 (mod 1). 

This simple classical picture is valid at low tempera- 
tures, ksT <C Eg, and high tunneling resistance, R t 3> 
Rk = h/e 2 ~ 25.8fc£l. Thermal fluctuations will round 
off the corners of the step structure. But even at zero 
temperature, electron tunneling processes can smear out 
the staircase, which for large tunneling approaches a 
straight line. These quantum fluctuations renormalize 
the ground state energy and lead to an effective charging 
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Several recent theoretical papers have addressed the 
renormalization of the charging energy due to electron 
tunneling processes. Perturbative expansions, valid for 
small values of the dimcnsionless tunneling conductance 
" = 5^*, Yield* 
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TABLE I: Analytical predictions for the leading correction to 
the exponential suppression, Eq. JHJ , of the effective charging 
energy in the limit of large a. 
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Ref. Method 
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^5 instanton 4ty a 

£ instanton 16ir 6 a :i 

perturbation theory 87r 6 « 3 (l — (5 log a)/(2ir 2 a)) 

,6 RG ~ a 

2 real-time RG 0.5(7ra) 6 ' 5 



In order to treat the opposite limit of large tunneling 
conductance, various approaches have been proposed - 
with controversial results. While these analytical predic- 
tions agree on the leading exponential suppression of the 
effective charging energy, 
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the pre-exponential factors differ widely. Using an instan- 
ton approach, Panyukov and Zaikin 3 obtain f(a) ~ a 2 , 
whereas Wang and Graberlj^ find f(a) ~ a 3 . Lukyanov 
and Zamolodchikov 5 use a perturbative expansion in 1/a 
and propose an expression similar to Ref. |j, with an ad- 
ditional logarithmic term. The renormalization group 
(RG) calculation by Hofstetter and Zwerges^ predicts a 
linear dependence, while the "real-time RG" approach of 
Konig and Schoelleui leads to a pre-factor f(a) ~ a 6 ' 5 
in the limit of large a. We list these different theoretical 
predictions in Tab. [I] 

Several attempts to resolve the issue by numerical 
simulations&ifliii have failed because of size constraints 
and the inability to reach the asymptotic regime where 
the exponential factor in Eq. © dominates. Some of 
these Monte Carlo results even seemed to disagree for 
intermediate values of a. 

The purpose of this study is to provide accurate data 
far into the strong tunneling regime, enabling us to test 
the predictions from the various analytical approaches. 
We observe two distinct behaviors, depending on the in- 
verse temperature (3 and the tunneling strength a. If 
(3Eq <c 1, multiple phase slip paths get suppressed 
and we find a pre-exponential factor which roughly 
agrees with the single-instanton predictions^ In the limit 
(3Eq 3> 1, which is relevant for the comparison with an- 
alytical results for zero temperature, we find f(a) ~ a 5 . 
This result disagrees with all the predictions in Tab. [i] 



II. MODEL AND METHOD 
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In Eq. (@J) , the sum is over all paths with winding number 
±n and the imaginary time effective action reads (H = 1) 

+ .//^^l (a, 

Jo Jo sin (|(r - t')) 

The first term accounts for the charging energy and the 
long ranged part for electron tunneling. 

This system can be mapped to a chain of classical XY 
spins, as the discretized action becomes 



j=i j<j' 
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We use a Trotter number N, corresponding to an imag- 
inary time step (and short-distance cutoff) At = /3/N 
and nearest neighbor coupling T = l/(2EcAr). Periodic 
boundary conditions <pN+i = 4>l are employed. 
In Fourier space the action becomes local, 
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where ipk — Ylj=i e * ™ jfe e 1 ^ denotes the Fourier trans- 
form of e 1 ^ and gk the Fourier transform of the kernel 

Cj + o) 

Using fast Fourier transformation 1 " 3 " it is possible to com- 
pute the action of a configuration in a time 0(N log N) 
despite the long ranged interactions. 

The winding number lu of a phase configuration 
{4>j}j=i,...N is determined by summing up the phase dif- 
ferences as 



JV 

w = ^2[(4>j+i - <t>j +tt) mod 2tt] - Nn. (9) 

The effective charging energy can then be computed from 
the expectation value of the winding number squared^, 



The partition function of the single electron box with 
gate charge zero can be written as a path integral over a 
compact angular variable cj> (conjugate to the number of 
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A. Cluster Monte Carlo 

In most of the following analysis, we use path in- 
tegral Monte Carlo simulations to study the action in 
Eq. |JBJ , employing a variant of the Swendsen-Wang clus- 
ter algorithmi 4 * 1 ^ to reduce auto-correlation times. The 
efficient treatment of long ranged interactions proposed 
by Luijten and Blote 1 ^ reduces the time for a cluster up- 
date to 0(N log N). This allows us to simulate chains of 
up to 10 7 spins, which is three orders of magnitude larger 
than the systems which have previously been studied. 

We ran simulations for AtEc = 0.5, 0.125, 0.05 and 
0.005 and found essentially the same results for the latter 
three values of the discretization step (especially for large 
a, which is the limit of interest). We therefore consider 
an extrapolation in At unnecessary for the purpose of 
this study and present data for AtEc = 0.05, unless 
otherwise stated. With this choice of At, the cluster 
algorithm allows us to simulate the single-electron box 
down to temperatures f3Ec ~ 10 4 . However, for large 
a, the fraction of paths with winding number different 
from zero decreases as exp(— it 2 a). Due to this rapidly 
deteriorating efficiency, the cluster algorithm can only be 
used in the range < a < 3. 



B. Transition Matrix Monte Carlo 

In order to compute the effective charging energy at 
even larger values of the tunneling conductance, we de- 
veloped a new Monte Carlo approach, which attempts to 
insert or remove phase slips (kinks) and in doing so mea- 
sures the relative weights of the different winding number 
sectors. As only paths of winding number and 1 con- 
tribute in the limit of large a, we restrict the discussion to 
the calculation of the relative weight of these two wind- 
ing number sectors. However, the algorithm can also be 
used to study additional winding number sectors. 

We denote by w n the weight of the configurations with 
winding number ±n, 

The effective charging energy (|1U|I can then be expressed 
for large a as 

Eh _ 2^ 2 E^ 2 ^ „ 2tt 2 m 

Ec (3E C En=o w n ~ PE c w ' 1 ' 

The weights wq and w\ correspond to a "density of 
states" in winding number space. In order to determine 
wi/uuq we could sample the winding number space with 
the inverse density of states and energy space with the 
usual Boltzmann weight. The acceptance probability for 



a proposed kink update from a configuration with wind- 
ing number n and action S to winding number n' and 
action S' then becomes 



P((n, S) -> On', S')) = min 1, — exp(S - S') , (13) 

V w n > J 

whereas a cluster update from n to n' would be accepted 
with probability 

P((n) -» fn')) = min ( 1, — ^1 . (14) 

V w n< J 

The relative density of states w\/wq then satisfies the 
"flat histogram" condition, which for n S {0, 1} can be 
expressed as 

(P(0 -> 1)} {wi/wo) = (P(l -> 0)} ( Wl /w ). (15) 

Averages are taken over a sequence of kink- and cluster 
updates. 

Our strategy is to determine the transition probabili- 
ties (P(0 — > 1)) and (P(l — > 0)) rather than trying to find 
the relative occupation of the winding sectors wi /wq by a 
Wang-Landau type iterative procedureit This is an ap- 
proach in the spirit of the transition Matrix Monte Carlo 
methodic. 

Each winding number sector is treated separately and 
the cluster updates serve essentially to randomize the 
configurations within that sector, although their contri- 
bution to the probabilities {P{i — > j)) becomes important 
for smaller a. In a kink update we insert a phase slip 

^ k A) =±2arctan(^) (16) 

of random width A at some random position To (or rather 
'/'kink'' = ±7rarctan((T — T )/A)/arctan(/3/2A), periodi- 
cally continued outside the the interval — <| < t — To < f , 
which is a slightly modified version compatible with the 
finite size of the system). These are stationary paths of 
the long ranged part of the action (JSJ) in the limit (3 — > oo. 

It is important to consider the entire range of widths 
including stretched out configurations. In order to study 
the typical distribution of widths we generated configu- 
rations for N — 10 5 ((3Ec — 5 • 10 3 ) in winding sector 
1 using the cluster algorithm and applied a decimation 
procedure in which the number of spins was reduced by 
half at each step. The average winding number (|a;( n )|) 
of these coarse-grained configurations is shown as a func- 
tion of decimation steps n in Fig. [31 Removing every 
second spin reduces the width of the kink to approxi- 
mately half its previous value until the kink disappears 
completely. From the dependence of (|a/ n )|) on n it is 
therefore possible to gain insight into the original dis- 
tribution of kink widths in winding sector 1. One finds 
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decimation step n 

FIG. 3: Average winding number (|ur n '|) of coarse grained 
configurations. The original configuration has winding num- 
ber = 1. It is renormalized by removing every second 
spin. Kinks of width « 2™ will disappear after approximately 
n iterations, resulting in a coarse-grained configuration of 
winding number 0. 



from Fig. |3 where we show the average winding after 
coarse-graining a configuration, that the kink widths for 
1.6 < a < 2.4 are approximately uniformly distributed in 
the range [2 6 ,2 16 ], whereas for a = 9 only smooth con- 
figurations with widths 6 [2 13 ,2 16 ] are generated. The 
sharpest kinks have a width of 2 5 . In order to cover the 
whole range of relevant paths, we chose A m i n = 10/3/iV 
and A max = 0. 

Inserting <%°J^ changes the original configuration <fi 
with winding number n to (j)' with winding number nil. 
The corresponding Boltzmann weights exp(S'[0] — ^I^']) 
are computed using Eq. J7J and stored during the simula- 
tion. Updates which produce a configuration with wind- 
ing number \n\ > 1 as well as cluster updates which leave 
the winding number unchanged get the weight 0. A clus- 
ter update which connects to the other winding number 
sector gets the weight 1. From these data, one can cal- 
culate the average transition probabilities (P(i — > j)} as 
a function of w\/wq using Eqs. 111311 and l|14fl . Finally, 
the relative weight of winding sector 1 is determined by 
solving Eq. l(ToT) . 

To illustrate this procedure, Fig. 0] shows the proba- 
bilities (P(0 -> 1)} (with positive slope) and (P(l -> 0)) 
(with negative slope) as a function of — \n.{w\ / 'wq) for sev- 
eral values of a. The intersection points of these curves 
determine wi/wq- The effective charging energies ob- 
tained by this new method are consistent with the results 
from cluster Monte Carlo for the values of a which can 
be treated with the latter method. The new approach, 
however, allows us to simulate the system at much higher 
tunneling conductance, as it improves the sensitivity to 
winding number fluctuations by dozens of orders of mag- 
nitude. 
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FIG. 4: Transition probabilities {P{0 — > 1)} (with positive 
slope) and (P(l — > 0)} (with negative slope) plotted as a func- 
tion of — ln(wi/tt)o). The intersection points of these curves 
determines wi/wo- From left to right, the data correspond to 
a = 2, 3, 9. 
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FIG. 5: Effective charging energy for f3Ec = 5 • 10 3 computed 
using cluster updates and kink updates. Previously published 
results obtained by means of local update schemes are shown 
for comparison. 



III. EFFECTIVE CHARGING ENERGY 

A. Leading exponential behavior 

The range of applicability of the different simulation 
methods is illustrated in Fig. [SI where we plot the effec- 
tive charging energy as a function of tunneling conduc- 
tance for j3Ec = 5 ■ 10 3 . While cluster updates allow to 
compute accurate data at small and intermediate values 
of a, they lose their effectiveness around Eq/Ec ~ 1CP 8 
or a ~ 3. On the other hand, the kink- updates enable us 
to cover tunneling strengths beyond a = 9, which corre- 
sponds to E* c jE c ~ 10~ 34 . This highly sensitive method 
allows us to clearly observe the theoretically predicted ex- 
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FIG. 6: Effective charging energy in the weak and inter- 
mediate tunneling regime. The black lines show the first-, 
second- and third-order perturbation results @. The data 
for (3Ec = 5 ■ 10 3 were computed with At = 0.005 in order 
to reduce the energy shift due to discretization effects. 
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FIG. 7: Correction to the exponential behavior, Eq. [3] at 
intermediate and large tunneling conductance and compari- 
son with analytical predictions. A crossover from intermedi- 
ate tunneling (zero temperature) behavior to strong tunneling 
(intermediate temperature) behavior occurs at (5Eq w 1. 



ponential suppression J2J of the effective charging energy. 
For comparison, we also plot some results from local up- 
date simulationsiiSiii Although these simulations were 
performed on much smaller lattices, they merely reached 
E* c /E c « lO" 2 . 



B. Weak and intermediate tunneling 

While the weak tunneling regime can be treated per- 
turbatively, few theoretical results exist for the experi- 
mentally accessible intermediate region. Since previous 
Monte Carlo simulationa&iSiii seemed to disagree in this 
region, we plot our data for weak and intermediate tun- 
neling in Fig. [S] The third order perturbation result @ 
reproduces our data for f3Ec — 5- 10 3 and At = 0.005 up 
to a = 0.6 (the data for At = 0.05 are shifted to slightly 
larger E* c jEc)- Also shown in Fig.0are results for lower 
and higher temperature and the data of Wang et ali& and 
Herrero et alM>. The results in Ref. ^] were obtained for 
f3Ec — 5 • 10 2 and are consistent with our data. The 
larger time step AtEq — 0.2 used in their simulation 
may be the reason for the tendency to somewhat larger 
values of E* c . The systematically lower values in Ref. fill 
were computed at higher temperature f3Ec = 10 2 , which 
explains the discrepancy. 

Around a « 0.6 the system enters an intermediate 
tunneling regime, where perturbation theory fails but 
the asymptotic formulas for large tunneling are not yet 
valid. The exponential suppression of the effective charg- 
ing emerges above a « 1.5 (see Figs. and EJ). This 
strong tunneling region cannot be accessed by local up- 
date schemes. 



C. Strong tunneling and crossover 

In order to determine the pre-exponential factor f(a) 
in Eq. @, we multiply the data with exp(7r 2 cv). The 
value 7r 2 a corresponds to the (long ranged) action of an 
optimal phase slip path l|16fl hi the limit (3 — > oo and 
At — > 0. Since the leading power of the pre-exponential 
factor will depend sensitively on this value, one might 
wonder whether the discretization leads to a significant 
change. For our values of (3 and At, however, we found 
that periodic versions of the paths (|16|) , with A not too 
small, have an action which deviates very little from 7r 2 a. 

In Fig. [7| we plot cluster Monte Carlo results for 
f3E c = 5 • 10 2 , 5 • 10 3 , 5 • 10 4 and 5 • 10 5 . The lowest 
temperature was computed with AtEc = 0.2, the oth- 
ers with AtEc = 0.05. Also shown in Fig. are the 
asymptotic predictions from several analytical calcula- 
tions. The line with slope 2 (in the log-log plot) is the 
result of Panyukov and Zaikir. 3 , those with slope 3 were 
predicted in Refs. and while the line with slope 6.5 
shows the asymptotic form of the "real time RG" result by 
Konig and Schoellei— . Hofstetter and Zwerger— predict a 
slope of 1, but not the amplitude. 

As the tunneling strength increases, the effective charg- 
ing energy Eq gets reduced and for fixed temperature one 
observes a crossover from intermediate tunneling behav- 
ior, where f3E* c 3> 1, to strong tunneling behavior in 
the region (3E^ <C 1, or equivalently for a fixed, from 
zero temperature behavior to an intermediate temper- 
ature behavior. This crossover occurs around (3E^ = 1 
and is clearly visible in Fig.0 For the temperatures plot- 
ted in this figure, the crossover condition corresponds to 
cv-values in the range 1.4 < a < 2.2. 
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FIG. 8: The pre-exponential factors in the strong tunneling 
regime can be fitted by shifted parabolas, Eq. 1181 . The ver- 
tical lines indicate the values of a, where the occupation of 
winding number sectors 1 and 2 vanish (see Fig. ©). Data 
points marked with circles show the results for /3Ec = 5 ■ 10 3 
computed by the transition matrix method, which considers 
configurations with kink-anti-kink pairs. 

D. Zero temperature regime, /3Eq 3> 1 

Theoretical predictions for zero temperature have to be 
compared to the Monte Carlo data in the zero temper- 
ature (intermediate tunneling) limit (3Ec ^> (3Eq 3> 1. 
Because of size restrictions we can only obtain accurate 
results for a < 1.6. These data points can be fitted with 
a power-law as indicated by the dotted line in Fig. [7| It 
corresponds to a pre-exponential factor 

f(a) ~a 5 . (17) 

Because this fit is performed in the region of interme- 
diate tunneling strengths, it is possible that the slope 
of the zero-temperature curve will decrease and eventu- 
ally approach for example the solution of Lukyanov and 
Zalomodchikov. - The calculation of Ref. Moverestimates 
the exponent, whereas the result of Ref. y lies too low. 
It is even more unlikely that the zero-temperature curve 
will approach the slope 1 predicted in Ref. @. The lat- 
ter RG-calculation was based on the assumption that the 
winding number remains unchanged under renormaliza- 
tion, which appears incorrect, as can be seen in Fig. [3] 

E. Intermediate temperature regime, /3Eq -C 1 

In the intermediate temperature (strong-tunneling) re- 
gion (3Ec > 1 > 0Eq, the pre-exponential factor can 
be fitted by a shifted parabola 

f{a) ~ a 2 + const. (18) 

In Fig.[S]we show these parabolic fits for (3Ec = 5 • 10 2 , 
5 ■ 10 3 and 5 • 10 4 . Upon close inspection, these data show 
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FIG. 9: Relative occupation of the winding number sectors 1, 
2 and 3 as a function of tunneling strength. As a increases, 
higher-winding-number configurations get suppressed and in 
a finite number of updates, multiple kink paths will not be 
generated at large values of a. 
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FIG. 10: Pre-exponential factor of the effective charging en- 
ergy in the intermediate and strong tunneling regime. The 
parabolic fit to the cluster Monte Carlo data (see Fig. |SJ is 
consistent with the data obtained from the transition matrix 
approach even at much higher tunneling conductance. 



some intriguing features. At certain values of a, the slope 
seems to change suddenly and the quality of the data 
deteriorates. This feature is a result of the suppression of 
multiple- kink paths. As a is increased, the occupation of 
higher winding sectors is reduced and for a finite number 
of measurements it will eventually vanish. The visible 
change in slope occurs at the value of a, where multiple- 
kink paths get completely suppressed. 

To illustrate this point, we plot in Fig.[5]the relative oc- 
cupation of the lowest three winding number sectors as a 
function of a. The position where the winding numbers 
2 and 3 get suppressed are marked by the vertical line 
segments in Figs. |H1 and |5J The transition matrix Monte 



7 



Carlo data for very large tunneling strength are consis- 
tent with the parabolic fit to the cluster Monte Carlo 
data, as shown in Fig. ED 

The instanton calculations^ neglect the interactions 
between phase-slips. This is probably the reason why 
they predict a pre-exponential factor which is roughly 
consistent with the simulation results in the large tun- 
neling (intermediat temperature) region, where multiple 
phase slips no longer contribute. In the zero temperature 
regime (3Eq 3> 1, however, multiple-kink configurations 
are abundant (see Fig. 0. The considerably different 
slope in this region leads to the conclusion, that inter- 
actions between phase slips are important. They should 
be considered in a theory which attempts to predict the 
pre-exponential factor f(a) for the single-electron box at 
zero temperature. 

IV. CONCLUSIONS 

We calculated the effective charging energy of a single- 
electron box using cluster Monte Carlo simulations and 
a new transition matrix Monte Carlo approach. The new 
method allows us to study very large tunneling strengths 
and to observe the leading exponential suppression of the 
charging energy, E* c jEc — f(a) exp(— ir 2 a), over more 



than thirty orders of magnitude. For a finite inverse tem- 
perature one observes two regimes. In the large tunnel- 
ing (intermediate temperature) regime, EqP <C 1, where 
multiple kink paths are suppressed, the pre-exponential 
factor has a leading power /(a) ~ a 2 , consistent with 
the instanton calculation of Ref. 0. In the intermediate 
tunneling (zero temperature) regime, E^(3 3> 1, mul- 
tiple kink paths are abundant and the interactions be- 
tween phase-slips become important. Simulation results 
obtained in this regime can be compared to analytical 
predictions for zero temperature. Fitting our data in 
the range 1.0 < a < 1.6 to a power-law, we find a pre- 
exponential factor f{a) ~ a 5 , which is not consistent 
with any of the theories. 
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